;;--------------------------------------------------------------------
;; Concepts illustrated:
;;   - Drawing box plots
;;   - Using a box plot to show the the median, minimum/maximum value, 
;;     and the 10th/90th percentiles of two timeseries
;;   - Adding text to a box plot
;;   - Sorting data
;;   - Setting the mininum/maximum value of the Y axis in a box plot
;;   - Generating dummy data using "random_normal"
;;--------------------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "$KAILIB_ROOT/callvar.ncl" 


begin

  nexp = 6 

  season = "DJF" 
  season = "SON" 
  season = "JJA" 
  season = "MAM" 
  siteid = 38   ;;capepoint
  passid = 9    ;;capegrim 
  H1215  = False 
  cpp    = "" 

  ss = (/ "ANN","DJF","JJA","MAM","SON" /)

  datapata = "$KAILIB_ROOT/sct_creatbox/mod/"  
  datapatb = "$KAILIB_ROOT/sct_creatbox/obs/"  
  datapatc = "$KAILIB_ROOT/sct_creatbox/"  

  fln_sitename = "radon_sitename.dat" 

  sname = asciiread(datapatc+fln_sitename,-1,"string")

  ns = dimsizes(sname)

;;-----------------------------------------------------------------
;; plot box plots for all sites avaiable 
;;-----------------------------------------------------------------
  

;;-----------------------------------------------------------------
;; load settings by user 
;;-----------------------------------------------------------------
  load "load_info.ncl" 

  print("") 
  print("season: "+season) 
  print("") 


  boxa = new((/ns,nexp,5/),"float") 
  boxb = new((/ns,nexp/),  "float") 
  boxc = new((/ns  /),  "float") 

;;-----------------------------------------------------------------
;; filename 
;;-----------------------------------------------------------------
  filenamb = datapata+"radon_mine_E0_1999-2003_"+season+".nc" 
  filenamc = datapata+"radon_mine_E1_1999-2003_"+season+".nc" 
  filenamd = datapata+"radon_mine_E2_1999-2003_"+season+".nc" 
  filename = datapata+"radon_mine_E8_1999-2003_"+season+".nc" 
  filenamf = datapata+"radon_mine_E0clim_1999-2003_"+season+".nc" 
  filenamg = datapata+"radon_mine_E3clim_1999-2003_"+season+".nc" 
  filenamx = datapatb+"radon_obs_mm.nc"
  filenamo = datapatc+"box_data_"+season+".nc" 

  system("rm "+filenamo) 

  flb = addfile(filenamb,"r") 
  flc = addfile(filenamc,"r") 
  fld = addfile(filenamd,"r") 
  fle = addfile(filename,"r") 
  flf = addfile(filenamf,"r") 
  flg = addfile(filenamg,"r") 
  flo = addfile(filenamo,"c") 

;;-----------------------------------------------------------------
;; loop over sites 
;;-----------------------------------------------------------------
  do is = 0,ns-1

  obssite = sname(is)  

  if(sname(is).eq."cape_grim") then 
     obssite = "capegrim" 
  end if 

  if(sname(is).eq."cape_point") then 
     obssite = "capepoint" 
  end if 

  filenama = datapatb+"radon_obs_"+obssite+"_"+season+".nc"

  haveobs = isfilepresent(filenama) 

  if(haveobs) then 
     print("open filename: "+filenama) 
     fla = addfile(filenama,"r")
  end if
  
  ;filenamm = obsfilepath+"radon.obs."+sitenama+".9999.dat"
  ;print(""+filenamm) 
  ;rmm = asciiread(filenamm,-1,"float") 
  ;print("plot without OBS data")

  flx = addfile(filenamx,"r")   
  rmmobs = flx->radon
  rmm = rmmobs(:,is) 
  rmm@_FillValue = -999
  rmma = (/1.,1.,1./) 

  if(season.eq."ANN") then 
     ram = dim_avg(rmm)
  end if 

  if(season.eq."DJF") then 
     rmma(0) = rmm(11) 
     rmma(1) = rmm(0) 
     rmma(2) = rmm(1) 
     ram = dim_avg(rmma) 
  end if 

  if(season.eq."JJA") then 
     rmma(0) = rmm(5) 
     rmma(1) = rmm(6) 
     rmma(2) = rmm(7) 
     ram = dim_avg(rmma) 
  end if 

  if(season.eq."MAM") then 
     rmma(0) = rmm(2) 
     rmma(1) = rmm(3) 
     rmma(2) = rmm(4) 
     ram = dim_avg(rmma) 
  end if 

  if(season.eq."SON") then 
     rmma(0) = rmm(8) 
     rmma(1) = rmm(9) 
     rmma(2) = rmm(10) 
     ram = dim_avg(rmma) 
  end if 


;;-----------------------------------------------------------------
;; read data 
;;-----------------------------------------------------------------
  xbb = flb->radon(:,is)  ;;rn_1m 
  xcc = flc->radon(:,is)  ;;rn_1m 
  xdd = fld->radon(:,is)  ;;rn_1m 
  xee = fle->radon(:,is)  ;;rn_1m 
  xff = flf->radon(:,is)  ;;rn_1m 
  xgg = flg->radon(:,is)  ;;rn_1m 

  if(haveobs) then 
     xab = fla->rn_1m 
  else 
     xab = xbb
     xab = 0. 
  end if 

;;-----------------------------------------------------------------
;; set negative values to missing 
;;-----------------------------------------------------------------
  xab = where(xab.gt.0.,xab,-999) 
  xbb = where(xbb.gt.0.,xbb,-999) 
  xcc = where(xcc.gt.0.,xcc,-999) 
  xdd = where(xdd.gt.0.,xdd,-999) 
  xee = where(xee.gt.0.,xee,-999) 
  xff = where(xff.gt.0.,xff,-999) 
  xgg = where(xgg.gt.0.,xgg,-999) 

;;-----------------------------------------------------------------
;; select obs data 
;;-----------------------------------------------------------------
  if(sname(is).eq."capegrim".and.H1215) then 
     xaa = var_1HH1215(xab,1)      ;; 1D to 2D 
     ybb = var_3HH1215(xbb,1)      ;; 1D to 2D 
     ycc = var_3HH1215(xcc,1)      ;; 1D to 2D 
     ydd = var_3HH1215(xdd,1)      ;; 1D to 2D 
     yee = var_3HH1215(xee,1)      ;; 1D to 2D 
     yff = var_3HH1215(xff,1)      ;; 1D to 2D 
     ygg = var_3HH1215(xgg,1)      ;; 1D to 2D 
     cpp = "H1215"
  else
     xaa = xab 
     ybb = xbb 
     ycc = xcc 
     ydd = xdd 
     yee = xee 
     yff = xff 
     ygg = xgg 
  end if 

  if(sname(is).eq."capegrim") then 
     passid = 9 
  end if
 
  if(haveobs) then 
  xaa = var_passend(xaa,passid)   ;; convert to GMT time  
  else
  xaa = 0.  
  end if
 
  yaa = xaa(::3)                  ;; only 3H data  

;;-----------------------------------------------------------------
;; remove missing values 
;;-----------------------------------------------------------------
  aa = vnmiss_1D(yaa) 
  bb = vnmiss_1D(ybb) 
  cc = vnmiss_1D(ycc) 
  dd = vnmiss_1D(ydd) 
  ee = vnmiss_1D(yee) 
  ff = vnmiss_1D(yff) 
  gg = vnmiss_1D(ygg) 



  dm = dimsizes(aa)

;;-----------------------------------------------------------------
;; get 10% 25% 50% 75% 90% percentile index 
;;-----------------------------------------------------------------
  x10 = round(.10*dm,3)-1  
  x25 = round(.25*dm,3)-1     ; -1 to account for NCL indexing starting
  x75 = round(.75*dm,3)-1     ; at 0            
  x90 = round(.90*dm,3)-1  

;;-----------------------------------------------------------------
;; sort the data 
;;-----------------------------------------------------------------
  qsort(aa)    
  qsort(bb)  
  qsort(cc)  
  qsort(dd)  
  qsort(ee)  
  qsort(ff)  
  qsort(gg)  

  pp=new((/nexp,5/),float) 
  pa=new((/nexp/),float) 
  pb=new((/1/),float)

  pp(0,:) = (/aa(x10),aa(x25),dim_median(aa),aa(x75),aa(x90)/)  

  if(haveobs) then 
  pa(0) = dim_avg(aa)  
  pb    = ram 
  else
  pa(0) = ram 
  pb    = ram 
  end if 

;;-----------------------------------------------------------------
;; do the same for model data 
;;-----------------------------------------------------------------
  dm = dimsizes(bb) 

  x10 = round(.10*dm,3)-1  
  x25 = round(.25*dm,3)-1  
  x75 = round(.75*dm,3)-1     ; at 0            
  x90 = round(.90*dm,3)-1  

  pp(1,:) = (/bb(x10),bb(x25),dim_median(bb),bb(x75),bb(x90)/)

  dm = dimsizes(dd) 

  x10 = round(.10*dm,3)-1  
  x25 = round(.25*dm,3)-1  
  x75 = round(.75*dm,3)-1     ; at 0            
  x90 = round(.90*dm,3)-1  

  pp(2,:) = (/dd(x10),dd(x25),dim_median(dd),dd(x75),dd(x90)/)

  dm = dimsizes(ee) 

  x10 = round(.10*dm,3)-1  
  x25 = round(.25*dm,3)-1  
  x75 = round(.75*dm,3)-1     ; at 0            
  x90 = round(.90*dm,3)-1  

  pp(3,:) = (/ee(x10),ee(x25),dim_median(ee),ee(x75),ee(x90)/)

  dm = dimsizes(ff)

  x10 = round(.10*dm,3)-1
  x25 = round(.25*dm,3)-1
  x75 = round(.75*dm,3)-1     ; at 0
  x90 = round(.90*dm,3)-1

  pp(4,:) = (/ff(x10),ff(x25),dim_median(ff),ff(x75),ff(x90)/)

  dm = dimsizes(gg)

  x10 = round(.10*dm,3)-1
  x25 = round(.25*dm,3)-1
  x75 = round(.75*dm,3)-1     ; at 0
  x90 = round(.90*dm,3)-1

  pp(5,:) = (/gg(x10),gg(x25),dim_median(gg),gg(x75),gg(x90)/)

  pa(1) = dim_avg(bb)  
  pa(2) = dim_avg(dd)  
  pa(3) = dim_avg(ee)  
  pa(4) = dim_avg(ff)  
  pa(5) = dim_avg(gg)  

  boxa(is,:,:) = pp(:,:) 
  boxb(is,:)   = pa(:) 
  boxc(is)     = pb 

;;-----------------------------------------------------------------
;; clean mem 
;;-----------------------------------------------------------------
  delete(aa) 
  delete(bb) 
  delete(cc) 
  delete(dd) 
  delete(ee) 
  delete(ff) 
  delete(gg) 
  delete(xaa) 
  delete(xab) 
  delete(xbb) 
  delete(xcc) 
  delete(xdd) 
  delete(xee) 
  delete(xff) 
  delete(xgg) 
  delete(yaa) 
  delete(ybb) 
  delete(ycc) 
  delete(ydd) 
  delete(yee) 

  end do ;; ns 
 
;;-----------------------------------------------------------------
;; defination 
;;-----------------------------------------------------------------
  boxa!0 = "site"
  boxa&site = ispan(1,ns,1) 
  boxa!1 = "exp"
  boxa&exp = ispan(1,nexp,1) 
  boxa!2 = "pt"
  boxa&pt = ispan(1,5,1) 

  boxb!0 = "site"
  boxb&site = ispan(1,ns,1)  
  boxb!1 = "exp"
  boxb&exp = ispan(1,nexp,1) 

  boxc!0 = "site"
  boxc&site = ispan(1,ns,1)  

;;-----------------------------------------------------------------
;; output 
;;-----------------------------------------------------------------
  flo->boxa=boxa
  flo->boxb=boxb
  flo->boxc=boxc

end



